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^\ • We implement the numerical method of summing Green function diagrams on the Matsubara 

0^ ' frequency axis for the fluctuation exchange (FLEX) approximation. Our method has previously been 

applied to the attractive Hubbard model for low density. Here we apply our numerical algorithm 
to the Hubbard model close to half filling (p — 0.40), and for T/t — 0.03, in order to study 
the dynamics of one- and two-particle Green functions. For the values of the chosen parameters 
we see the formation of three branches which we associate with the a two-peak structure in the 
imaginary part of the self-energy. From the imaginary part of the self-energy we conclude that our 
system is a Fermi liquid (for the temperature investigated here), since ImS{k,u}) « w around the 
chemical potential. We have compared our fully self-consistent FLEX solutions with a lower order 
approximation where the internal Green functions are approximated by free Green functions. These 
two approches, i.e., the fully selfconsistent and the non-selfconsistent ones give different results for 
^-v ' the parameters considered here. However, they have similar global results for small densities. 

^ ■ Pacs numbers: 74.20.-Fg, 74.10.-z,74.60.-w, 74.72.-h 
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High-temperature superconductors a wide range of behaviour atypical |^ of the standard band-theory of metals, 
since at half filling they should be metals while they happen to be insulators. Then, correlations are important. This 
Q ' leads us to consider that the theory of strongly correlated Fermin systems plays an important role for the pairing 
Q mechanism and other properties of the cuprates [|| . The current strategy (due to Anderson) to address the problem 
^ of high-Tc superconductivity is to try to find a theory accounting for the normal state properties of the cuprates. This 
• I— I , goes in analogy with the ordinary BCS theory where firstly the normal state is identified to be a Landau Fermi liquid 
' metal and then is found an electron pairing mechanism (unknown up to now in the cuprates) which destabilizes the 
normal state phase towards a superconducting state. In particular, due to their pronounced two-dimensional character 
- - the one-band Hubbard Hamiltonian of the Cu-3d hole states has been taken as one of the essential models |^ . The 
Hubbard model is the simplest many-particle model one can write down, which cannot be reduced to a single-particle 
theory Q. To study this model, the simultaneous evaluation of the frequency and momentum dependence of the 
one- and two-particle Green functions is crucial. Here, we adopt the opposite view of infinite dimensions which 
neglects the important role of spatial correlations. 

The Hubbard model has been a subject of intense numerical studies. In spite of exact diagonalization and quantum 
Monte Carlo studies ]9|-|l^], these methods are restricted to indeed small finite size systems. Another route which 
takes both dynamical as well as momentum space correlation into account is the use of diagrammatic approaches 
We should mention the self-consistent summation of all bubble and ladder diagrams, so called fiuctuation 



exchange approximation, or FLEX [^-18|, ||T^. This approximation (FLEX) is conserving in the sense of Baym and 
Kadanoff pl[|, i.e., it is consistent with microscopic conservation laws for particle number, energy, and momentum. 
We recall that the FLEX approximation belongs to the class of ^-derivable approximations which have been discussed 
a long time ago by Wortis in the context of perturbative series expansion of thermodynamic properties of models 
described by lattice Hamiltonians, mainly the Ising model. Wortis gives a set of steps to construct a ^-derivable 
approximation. 

In the FLEX approximation there is an interesting feature in the effective interaction: it is temperature and 
doping dependent through the spin susceptibility, x{Qt^)- The effective interaction depends on the properties of the 
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quasiparticles. One thus has a system in which, since the effective interaction both modifies the quasi-particle behavior 
and is itself alterated by that changed quasi-particle behavior, non-linear feedback, either positive or negative, can 
play a significant role |^^. According to the FLEX approach, the dominant contribution to the magnetic interaction 
between planar quasi-particles is assumed to come from spin-fluctuation exchange, and so will be proportional to 
x{q,uj). In the scheme of Pines and co-workers pBf , x('f, ^) is taken from experiments, while in the present work we 
have a fully self-consistent close set of equations to solve. 

Recently, Schmalian, Langer, Grabowski and Bennemann p^ ] have performed the calculation of the dynamical 
properties on the real axis, instead of the Matsubara frequencies. Here we perform our calculations on the Matsubara 
axis (1024 frequencies) obtaining the dynamical properties for all frequency range after analytically continuing the 
one- and two-particle Green functions with Pade approximants [ p5| . 

Our method for solving numerically the FLEX equations is based on the Fast Fourier tranformation (EFT) p^ ] 
which is suitable for the evaluation of the corresponding integrals. We refer to reference | |2^ ] where this technique 
is documented. In order to reach high densities, while working on the imaginary time, we have started from small 
densities (or low chemical potential). For each value of the chemical potential we have reached self-consistency and 
then we move in small steps to a higher chemical potential. We have proceeded this way since the next Green function 
is constructed from the nearby previous one. So, we always stay close to the real solution until we reach the density 
we are looking for. 

In Section O we describe our model Hamiltonian and the FLEX equations to be solved on the imaginary time. In 
Section [II, we present the data for the one- and two-particle spectral functions, the self-energy and density of states. 
It is observed the formation of three branches in the spectral density. We analyze, following the self-consistency of our 
equations, the presence of these two peaks as due to the presence of two peaks in the imaginary part of the effective 
interaction. We have compared our solutions with a low order approximation which consists in approximating the 
internal Green functions by free ones. Our results show that the calculations agree globally with each other for small 
densities. The dynamical quantities are radically different for the density/spin investigated in this paper (p = 0.40). 



In Section IV we discuss our result and present our conclusions. 



II. THE MODEL AND THE FLEX EQUATIONS. 

The Hubbard Hamiltonian is defined as 

H = -t ^ cj^ci'a + U^n^nii - fi^nia- , (1) 

<ll'>(j I la 

where the c\^{cia) are creation (annihilation) operators for electrons with spin a. The number-operator is nj^. = cj^cicr, 
t is the hopping matrix element between nearest neighbours I and I' , U is the on-site interaction and fi is the chemical 
potential in the grand canonical ensemble. Here we consider a repulsive interaction, U > 0. 
The one-particle Green function is expressed in terms of the self-energy, I](fc,iw„), by 

G{k,zcu„) = - ^ , (2) 

where e-^ = —2t{cos{kx) + cos{ky)) — /i. The self-energy of the Hubbard model (Eq. (|l|)) within the FLEX approxi- 
mation is given as follows |l4[ 

i;(fc,iw„) = — ^ yp_/i(q,ie„,.)G(k - q, iLj„ - ie^) , (3) 

where the effective interaction, V^-h(q, *em), resulting from the summation of the bubble and ladder diagrams (particle 
- hole channel), is given by 



Vp-hW, icm) = — x(q, iem) -j m : — T + 1 I rr / : T " 2 ) . (4) 



Here, 

T 



X(q, ie,n) ^ 'Yl + '5' 'i-<^m)G{k, iuj„) (5) 



k 
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is the particle-hole bubble with the renormalized Green function G(fc, itOn)- 'jJn = (2n + l)7rT and e„j = 2rmTT are the 
fermionic (odd) and bosonic (even) Matsubara frequencies. Let's state that the bubble and ladder sums represent the 
charge and longitudinal and the transversal spin fluctuations, respectively. The effect of particle-particle fluctuations 
(Cooper channel) has been left out of the present study. 

The previous set of equations closes with the expression for the p given by 

PmA^) = ^ + ^5]G(fczc^„) . (6) 

The form of Eq. (^ is well suited for numerical calculations. It has been justified by Schafroth, Rodrfguez-Nunez 
and Beck Eq. (^ represents a generalization of Mahan's Eq. (3.1.2) Eqs. (|[^) constitute the standard set 
of FLEX equations which have to be solved selfconsistently. 

In order to obtain results which are independent of finite size, one should use at least some 10^ Matsubara frequencies 
and a grid of 30 x 30 lattice points. The above scheme works in principle but a closer look at the equations for x and 
S shows that the straight forward implementation of these equations does not work in practice. This is due to the 
4~-fold loops which would occur in the computer program. Suppose we use 2000 Matsubara frequencies and a 30 x 30 
grid. Then we have to carry out for every grid point and every frequency the double sum over all frequencies and all 
grid points. This are of the order (30^ x 2000)^ = 3.24 x 10^^ complex operations. Even a typical super-computer 
would need one to several hours to make one iteration step. 

Since the frequency and momentum summations are convolutions, we evaluate them using the Fast Fourier Trans- 
form (FFT). In order to do that we need to transform the FLEX equations to direct space-time. For example, Eqs. 
(^,||) adopt the following forms 

x(x,T) = -G(x,r)G(x,-r) , (7) 

and 

S(x,r) = T/p_,,(x,r)G(x,T) . (8) 

In order be self-consistent, we calculate the frequency integrals on the imaginary time following a different approch 
to the one of Schmalian et al p^ . According to our experience, it is much better to work with Matsubara frequencies 
due to the fact that we have to handle neither Fermi nor Bose distribution functions. As we have previously said 
in the Introduction we have reached high densities by moving from a selfconsistent solution to another close by. In 



other words, we do not need going to real frequencies to have stability in our algorithm. In Section [Qig we present our 
results, making special emphasis on the the way we have reached high densities (see Fig. 1). 



III. NUMERICAL RESULTS AND INTERPRETATION 



In Fig. 1 we present the density/spin, p, as function of chemical potential, /i, for U/t — 4.0 and a temperature 
of T /t = 0.03. (This temperature is slightly above the temperature used in Ref. |24| of Schmalian et al. We have 
chosen this temperature since we can recover the asymptotic behavior, i.e., tt * 512 *T = A5t w 11.1J7. For T/t < 0.03 
we need to increase the number of Matsubara frequencies (we would need to go to RAM memories of 100 Mb or 
higher) in order to reach the asymptotic behavior. This asymptotic behavior is needed in order to compute the Fast 
Fourier transforms. In this figure we also show the the equivalent results for a low order approximation which consists 
in replacing in Eqs. (§-|^) the full one-particle Green function, G(fc, iw„), by the one-particle free Green function, 
Go{k,iuJn) — l/{iujn — £^)- This approximation was proposed by Serene and Hess |p9[| . We observe that the two 
curves almost coincide for low densities. For high densities, this is not the case, as we will see shortly. At this 
point we would like to emphasize that the way we have reached high densities settles the apparent drawbacks of the 
calculations of the T-matrix approach which is a lower order approximation to FLEX. There is recent calculation by 
Kagan et al where they have performed analytical calculations for the negative Hubbard model in the frame of 
the T— Matrix approximation. 

In Fig. 2 we present the spectral density, A^kjUj), as function of frequency, to for different momenta along the 
diagonal of the Brillouin zone (k, k) = 7r/16(n, n), with n an integer. We have the parameters U/t = 4.0, T/t = 0.03, 
p = 0.40. We are working with a system of 32 x 32 x 1024. The function A{k,uj) is obtained from the one-particle 
Green function by means of the analytical continuation, 

A(k,oj) = -- lim Im\G(k,uj + iS)] . (9) 

TT 5-^+0 
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In Eq. (||), we have chosen 5/t = 0.01. After selfconsistency, the chemical potential turns out to be ^ = — 0.95. 
We should like to point out that the chemical potential is independent of 6. This value is used at the end of the 
calculations to perform the analytical continuation as done by Vidberg and Serene See also Ref. |24| for more 
details. From Fig. 2 we see that the spectral functions have a well defined peak for small momenta and a satellite for 
positive frequency. For large momenta, we observe a kind of incoherent peak for positive frequencies and a satellite for 
negative frequencies. This incoherent peak could be decomposed into two peaks, one for small positive frequencies and 
another one for large positive frequencies. These results are different with respect to the negative Hubbard model at 
small density (T-Matrix Approximation), since in the T-Matrix approximation we find two peaks in the one-particle 
spectral function, while here we are resolving three peaks, instead. The relevance of the spectral functions in the 
Hubbard model both above and below the critical temperature, Tc of the superconducting state has been treated 
by Dahm and Tewordt (See their Figs. 17). They have compared their spectra with photoemission data of 
La2-xSrxCuO'i. 

For comparison, in Fig. 3, we present the spectral densities for the low order approximation. We observe that the 
peak structure is more pronunced than in the interacting case. This effect is also seen in the self-energy which we 
show next. 

In Fig. 4 we display the imaginary part of the self-energy, — /m[S(n, n, w)], along the diagonal of the Brillouin 
zone. The parameters are the same as in Fig. 2. We observe that energy dependence of — /m[S(n, n, w)], close to 
the chemical potential (the chemical potential is located at w = 0), has a gap-like structure. — /TO[S(n, n, w)] has a 
well defined gap structure close to the chemical potential (w « 0). We also see two almost symmetric peaks for every 
moment (fc, k). This presents a difference with respect to the attractive Hubbard model at small densities. For the 
attractive Hubbard model at small densities, there was a single peak in frequency for all values of (fc, k). This single 
peak structure in — /m[S(n, n, w)] was fitted by a two-pole Ansatz for the one-particle Green function. From Fig. 4 
we conclude that, for the studied temperature, we have a Fermi liquid system, since Iml](fc,cj) k, llP' for lo close to 
the chemical potential. We agree with Wermbtex |Q even when his calculation is mean-field like. See the discussion 
of Fig. 5. 

In Fig. 5 (a) we present the same quantity, — /TO[I](n, n, w)] vs u>, for the low order approximation. We inmediately 
see a well defined gap for high momenta. In consequence, the effect of fiuctuations is to reduce the exaggerated gap 
of the non-selfconsistent solution. Replacing our internal lines by free Green functions is equivalent to stay close to 
mean- field treatments, as it has been done by Kampf and Schrieffer pT| . These authors find well defined gaps in the 
spectral density which are most likely due to their mean field treatment. The approach mentioned in Ref. ||3l|| has 
been critically studied by Monthoux In Fig. 5 (b) we show — /to[S(7i, 0, cj)] vs uj. Here the gap is present for all 
the momenta. 

As it has been said in the Introduction the effective potential is strongly momentum and frequency dependent. We 
show in Figs. 6 and 7 the imaginary and real parts of the effective potential, respectively, for the same set of parameters 
as given in Fig. 2 along the diagonal of the Brillouin zone. We observe that Im\Vp-h{m,m,ijj)\ {Re\Vp-h{m,m,uj)\) 
is basically odd (even) in frequency around the chemical potential. This symmetry around zero is due to the fact that 
the Hartree shift (HS) has been substracted from the effective interaction. The real part of the effective interaction, 
i.e., i?e[V^_/i(m, 0, w)] vs lo compares qualitately well with Fig. 19(a) of Ref. The authors of Ref. (l9) have 

plotted the irreducible spin susceptibility, xi^i^) vs lo for some values of k. 

Now, we are in a position to explain the double peak structure in the imaginary part of the self-energy. Going to 
real frequencies, Eq. (|^) can be rewritten as 

if /■+°° 

i;(fc,z) = — ^ / dujVp-h{q.,z ~ uj)A{k - q,uj)nF{i^) 
q 1-"'-°° 

— Im[Vp^h[q,uj)]nB{Lo)G{q - k,-uj + z) , (10) 

where npiuj) and nsiLo) are the Fermi and Bose distribution functions, respectively. 

As the peaks of Im[Vp-h{q,to)] are symmetric in energies, then both terms of Eq. ( p^ ) contribute. As Im[Vp-h{q,Lo)] 
is antisymmetrical for frequencies close to /x where A{k, lo) has a peak, then we have two contributions to Im[T,{k, lo)]. 
Indeed, the numerical results for any (fc, k) of the self-energy mirror the behavior of Vp^hiQi ^)- Thus, the two-particle 
spectrum is introduced into the one-particle quantities like G{k,Lo) and A{k,uj). 
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IV. CONCLUSIONS. 



The FLEX approximation, a Baym-Kadanoff generalization of Hartree-Fock theory, has been implemented to study 
the frequency and momentum dependence of the one- and two-particle correlation functions. We have found the 
existence of three peaks in the one-particle spectral density, for the set of parameters investigated here. The presence 
of this structure in the spectral functions is a clear manifestation that correlations are indeed important. For smaller 
density/spin, for example p = 0.1 (Fig. 8), we find that the spectral functions are almost single peaks, or free-like 
quasi-particles. We have compared our self-consistent calculation with a low order approximation which consists in 
replacing the internal one-particle Green functions by free ones. The effects of the latter approximation are evident: 
the spectral functions are almost delta functions and the self-energy shows a wide gap around the chemical potential, 
which signals that this approximation is mean-field-like. We add that besides the good features present in the FLEX 
approximation [Q, we have shed some light on another aspect of it, i.e., the dependence of the one-particle properties 
on the two-particle Green functions due to the self-consistency of the FLEX equations. We mention the work of 
Nakamura, Moriya and Ueda and references therein where they point out the role of both the low and high 
frequency behavior of the spin fluctuations in the superconducting phase. For densities around half-filling correlations 
start to build in giving rise to three peaks in the energy spectrum. We argue that this analysis can be described in 
a generalized scheme of three-pole Ansatz for the spectral function |^ . Work along these line is in progress |36j in 
order to really resolve the three peaks in the one-particle spectral function. We should add that three peaks in the 
A{k,Lij) is equivalent to two peaks in the imaginary part of the self-energy. In this work, we have also discussed the 
differences of FLEX with respect to the T-Matrix approximation. 
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Figures. 

FIG. 1. The density/spin, p as function of chemical potential, p for the full self-consistent FLEXC equations as given in Eqs. 
(^,^). We also show the results for a low order approximation where we substitute the internal one-particle Green functions by 
free Green functions. The parameters used here are U/t = 4.0, T/t = 0.03. 

FIG. 2. The diagonal one-particle spectral function, A{n{TT , n) , lo) vs uj for different momenta along the diagonal of the 
Brillouin zone {k — (n,n)7r/16) for U/t = 4.0, T/t = 0.03. We have used an external damping of 5/t — 0.01, 32 x 32 points 
in the Brillouin zone and 1024 Matsubara frequencies. After self-consistent calculation of the coupled non-linear equations, we 
get p/t — —0.9558. We have runned our source code in single precision requiring 43 AIB of RAM memory. Each iteration 
takes 3.2 minutes of CPU time in a Pentium 166. 

FIG. 3. The diagonal one-particle spectral function, A{n{TT,n),ui) vs ui for different momenta along the diagonal of the 
Brillouin zone {k — (n, n)7r/16) for the low order approach. Same parameters as in Fig. 2. Here we have taken S/t = 0.0001. 

FIG. 4. — 7m[E(n(7r/16, 7r/16), tj)] vs uj for different momenta along the diagonal of the Brillouin zone (fc — {n,n)n/16). 
Same parameters as in Fig. 2. 
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FIG. 5. (a) — 7m[S(n(7r/16, 7r/16), w)] vs u) for different momenta along the diagonal of the Brillouin zone (fe = (n, n)7r/16) 
for the low order approximation. Same parameters as in Fig. 2. 5/1 = 0.0001. (b) — 7m[S(n(7r/16, o, w)] with the same 
parameters. 

FIG. 6. rm\Vp-h{jn{Ti /lQ),u:)] vs to for different momenta along the diagonal of the Brillouin zone [q — (m, m)7r/16). 
Same parameters as in Fig. 2. To simplify a little bit the notation, we have identified the effective potential, 
14)_/j(m(7r/16, 7r/16), w)] by T(m(7r/16, 7r/16),a;)]. The same is done in the next figure. 

FIG. 7. -Re[Vp-h(m(7r/16,7r/16),a')] vs u! for different momenta along the diagonal of the Brillouin zone (q = (m,m)7r/16). 
Same parameters as in Fig. 2. 

FIG. 8. The diagonal one-particle spectral function, A{n{TT , n) , u>) vs ui for different momenta along the diagonal of the 
Brillouin zone (k = {n, n)7r/16). Here the density/spin is 0.1 and n = —3.17t. 
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